Agent Based Models in Mata
Maarten Buis
University of Konstanz
maarten.buis@uni.kn
-------------------------------------------------------------------------------
index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- What is an Agent Based Model?
-------------------------------------------------------------------------------
What is an Agent Based Model?
An Agent Based Model is a simulation in which agents, that each follow
simple rules, interact with one another and thus produce a often
surprising outcome at the macro level.
The purpose of an ABM is to explore mechanisms through which actions of
the individual agents add up to a macro outcome, by varying the rules
that agents have to follow or varying with whom the agent can interact.
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- What is an Agent Based Model?
-------------------------------------------------------------------------------
Example 1: Immunization rate
Below is an example Agent Based Model that illustrates how a sufficiently
high immunization rate protects even the non-immunized individuals in a
population, and also allows one to explore at what immunization rate this
protection disappears.
In this model a single agent lives in each cell of a square grid and does
not move. An agent could have 5 states. She could be:
immunized
not-immunized and vulnerable
not-immunized and infected
not-immunized and infectious
not-immunized, healed, and no longer vulnerable
In the first round a single random vulnerable agent becomes infected.
In all subsequent rounds
all infected agents become infectious, and
infect all vulnerable agents within radius cell radius.
All infectious agents become healed and no longer vulnerable.
. set seed 123456789
. set rmsg on
r; t=0.00 16:52:06
. clear mata
r; t=0.00 16:52:06
. drop _all
r; t=0.00 16:52:06
. run infection.mata
r; t=0.01 16:52:06
.
. // run a single model and graph the evolution of the infection over time
. mata
------------------------------------------------- mata (type end to exit) -----
: foo = immune()
: foo.rdim(40)
: foo.cdim(40)
: foo.tdim(24)
: foo.radius(2)
: foo.not_vaccinated(240)
: foo.run()
: foo.extract(("id", "time", "status", "y", "x"))
: end
-------------------------------------------------------------------------------
r; t=0.15 16:52:06
.
. replace y = - y
(6,000 real changes made)
r; t=0.00 16:52:06
.
. mata st_local("rdim", strofreal(foo.rdim()))
r; t=0.00 16:52:06
. mata st_local("cdim", strofreal(foo.cdim()))
r; t=0.00 16:52:06
. mata st_local("tdim", strofreal(foo.tdim()))
r; t=0.00 16:52:06
. tempfile data
r; t=0.00 16:52:06
. save `data'
file C:\Users\Admin\AppData\Local\Temp\ST_2310_000008.tmp saved
r; t=0.00 16:52:06
. drop _all
r; t=0.00 16:52:06
. set obs `rdim'
number of observations (_N) was 0, now 40
r; t=0.00 16:52:06
. gen x = _n
r; t=0.00 16:52:06
. expand `cdim'
(1,560 observations created)
r; t=0.00 16:52:06
. bys x : gen y = -_n
r; t=0.00 16:52:06
. expand `tdim' + 1
(38,400 observations created)
r; t=0.00 16:52:06
. bys y x : gen time = _n
r; t=0.06 16:52:07
. merge 1:1 x y time using `data'
Result # of obs.
-----------------------------------------
not matched 34,000
from master 34,000 (_merge==1)
from using 0 (_merge==2)
matched 6,000 (_merge==3)
-----------------------------------------
r; t=0.08 16:52:07
. assert _merge != 2
r; t=0.00 16:52:07
. drop _merge
r; t=0.00 16:52:07
.
. forvalues t = 1/25 {
2. scatter y x if status == 0 & time == `t' , ///
> msymbol(Oh) || ///
> scatter y x if status == 1 & time == `t' , ///
> msymbol(S) mcolor(gs8) || ///
> scatter y x if status == 2 & time == `t' , ///
> msymbol(S) || ///
> scatter y x if status == 3 & time == `t' , ///
> msymbol(+) || ///
> scatter y x if status == . & time == `t', ///
> msymbol(o) msize(*.10) mcolor(gs8) ///
> aspect(1) name(gr`t', replace) ///
> yscale(off range(-41.5 -0.5)) ///
> xscale(off range(.5 41.5)) ///
> plotregion(margin(zero)) ///
> title("Iteration `t'") ///
> legend(order( 1 "vulnerable" ///
> 2 "infected" ///
> 3 "infectious" ///
> 4 "healed/immunized" ///
> 5 "immunized"))
3. }
r; t=16.13 16:52:23
.
. // run multiple models and plot the seriousness of the infection against
. // the immunization rate
. mata
------------------------------------------------- mata (type end to exit) -----
: scen = (16, 80, 120, 160, 200, 240, 280, 320, 360, 400)
: res = J(200,2,.)
: k=0
: for (i = 1; i <= 10 ; i++) {
> printf("number of vulnerable: %3.0f \n", scen[i])
> foo.not_vaccinated(scen[i])
> for(j = 1; j<= 20; j++) {
> k = k +1
> foo.run()
> res[k,.] = scen[i], foo.n_infected()
> }
> }
number of vulnerable: 16
number of vulnerable: 80
number of vulnerable: 120
number of vulnerable: 160
number of vulnerable: 200
number of vulnerable: 240
number of vulnerable: 280
number of vulnerable: 320
number of vulnerable: 360
number of vulnerable: 400
:
: idx = st_addvar("float", ("N","n"))
: st_addobs(200)
: st_store((1..200)',idx,res)
: end
-------------------------------------------------------------------------------
r; t=17.44 16:52:40
.
. gen perc_notvac = N/(40*40)*100
(40,000 missing values generated)
r; t=0.00 16:52:40
. gen perc_inf = n/N*100
(40,000 missing values generated)
r; t=0.00 16:52:40
.
. scatter perc_inf perc_notvac, msymbol(o) mcolor(%20) ///
> xtitle("percentage not vaccinated") ///
> ytitle("percentage infected of not vacinated ") ///
> name(scenarios, replace)
r; t=0.36 16:52:41
.
.
-------------------------------------------------------------------------------
Introduction -- What is an Agent Based Model?
-------------------------------------------------------------------------------
Example 2: Schellings segregation model
In this model there are blue and red agents, that live in cells on a
square grid. Each cell can contain only one agent, but not all cells are
occupied. An agent is happy if the proportion of agents of the same
colour in her neighbourhood is above a certain cut-off value. If she is
unhappy she will move to the nearest available cell that does satisfy her
needs.
The point of this simulation that even with a low cut-off value, like
30%, you will get fairly seggregated neighbourhoods.
. run schelling.mata
.
. mata
------------------------------------------------- mata (type end to exit) -----
: foo = schelling()
: foo.rdim(10)
: foo.cdim(10)
: foo.tdim(5)
: foo.nred(45)
: foo.nblue(45)
: foo.limit(.3)
: foo.run()
: foo.summtable()
Average proportion of neighbours with same color and average proportion of
happy agents
-----------------------------------
Iteration prop same prop happy
-----------------------------------
0 .556 .822
1 .758 .956
2 .799 .989
3 .808 .989
4 .82 .989
5 .829 1
-----------------------------------
: foo.extract(("r","c","t", "id", "red", "happy"))
: end
-------------------------------------------------------------------------------
.
. replace r = -r
(540 real changes made)
.
. forvalues t = 0/5 {
2. twoway scatter r c if red == 0 & t == `t' & happy == 1, ///
> msymbol(O) mcolor(blue) msize(*2) ///
> mlabel(id) mlabpos(0) || ///
> scatter r c if red == 1 & t == `t' & happy == 1, ///
> msymbol(O) mcolor(red) msize(*2) ///
> mlabel(id) mlabpos(0) || ///
> scatter r c if red == 0 & t == `t' & happy == 0, ///
> msymbol(Oh) mcolor(blue) msize(*2) ///
> mlabel(id) mlabpos(0) || ///
> scatter r c if red == 1 & t == `t' & happy == 0, ///
> title(iteration `i') ///
> msymbol(Oh) mcolor(red) msize(*2) ///
> mlabel(id) mlabpos(0) ///
> yscale(off range(-10.5 -.5)) ///
> xscale(off range(.5 10.5)) ///
> ylab(none) xlab(none) ///
> plotregion(margin(zero)) ///
> aspect(1) ///
> xline(.5(1)9.5) ///
> yline(-.5(-1)-9.5) ///
> legend(off) ///
> name(grsch`t', replace) title(Iteration `t')
3. }
-------------------------------------------------------------------------------
Introduction -- What is an Agent Based Model?
-------------------------------------------------------------------------------
Example 3: Bots influencing the mood
In this simulation there are agents that have a mood ranging from 10
(good) to 1 (bad).
In each iteration an agent reads a single newsitems. If the newsitem is
good it increases the mood, and if it is bad it decreases the mood.
The chance of choosing an item is proportional to the number of agents
that have previously read that article.
An item disappears after a certain number of iterations and a new
newsitem is created with a random number of initial "clicks".
A "bot" tries to decrease the mood by adding additional clicks to the
initial clicks of bad newsitems.
. run abm_clicks.mata
.
. mata
------------------------------------------------- mata (type end to exit) -----
: foo = clicks()
: foo.tdim(200)
: foo.pr_good(.5)
: foo.botclicks(50)
: foo.run()
: foo.mood(("m1","m2","m3","m4","m5","m6","m7","m8","m9","m10"))
: end
-------------------------------------------------------------------------------
. gen t = _n-1
. sort t
. forvalues i = 2/10 {
2. local j = `i' - 1
3. replace m`i' = m`i' + m`j'
4. }
(201 real changes made)
(201 real changes made)
(201 real changes made)
(201 real changes made)
(201 real changes made)
(201 real changes made)
(201 real changes made)
(201 real changes made)
(201 real changes made)
.
. twoway area m10 t, color("165 0 38") || ///
> area m9 t, color("215 48 39") || ///
> area m8 t, color("244 109 67") || ///
> area m7 t, color("253 174 97") || ///
> area m6 t, color("254 224 144") || ///
> area m5 t, color("224 243 248") || ///
> area m4 t, color("171 217 233") || ///
> area m3 t, color("116 173 209") || ///
> area m2 t, color(" 69 117 180") || ///
> area m1 t, color(" 49 54 149") ///
> plotregion(margin(zero)) ///
> name(mood, replace) ///
> legend(order(1 "10 good" ///
> 2 "9" ///
> 3 "8" ///
> 4 "7" ///
> 5 "6" ///
> 6 "5" ///
> 7 "4" ///
> 8 "3" ///
> 9 "2" ///
> 10 "1 bad" ))
. drop _all
. mata foo.items(("id", "t", "what", "val"))
. reshape wide val, i(id t) j(what)
(note: j = 1 2)
Data long -> wide
-----------------------------------------------------------------------------
Number of obs. 8040 -> 4020
Number of variables 4 -> 4
j variable (2 values) what -> (dropped)
xij variables:
val -> val1 val2
-----------------------------------------------------------------------------
.
. gen clicks_good = val2*val1
. gen clicks_bad = val2*(1-val1)
. collapse (sum) clicks_good clicks_bad, by(t)
. twoway line clicks_good clicks_bad t, ///
> lcolor("165 0 38" " 49 54 149") ///
> lpattern(solid solid) ///
> legend(order(1 "good" 2 "bad")) ///
> ytitle(total number of clicks) ///
> name(clicks, replace)
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- What is an Agent Based Model?
-------------------------------------------------------------------------------
What is this talk about?
This talk reports on an experiment to find out if it is viable to do ABMs
in Stata.
In particular many ABMs have tasks in common, e.g. in many cases agents
"live" inside a square grid and we want to know who their "neighbours"
are. How can we create and maintain a library of these common functions.
This library needs to be:
convenient to "import" in your model
modular
easy to maintain and extend
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- What is an Agent Based Model?
-------------------------------------------------------------------------------
The alternatives
There are various dedicated environments for creating ABMs.
netlogo (netlogo) http://ccl.northwestern.edu/netlogo/
repast (java, python, C++, and others) https://repast.github.io/
mason (java) https://cs.gmu.edu/~eclab/projects/mason/
mesa (python) https://mesa.readthedocs.io/
People who primarily create ABMs are not going to move away from those,
and rightly so.
So the toolkits discussed in this presentation should be aimed at people
who primarily use Stata and occationally want to make an ABM.
The question becomes: can we make making ABMs in Stata
easier/quicker/more convenient for this target audience than learning one
of these environments.
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- Mata
-------------------------------------------------------------------------------
Entering and exiting Mata, scalars, matrices
For implementing ABMs we will be using Stata's matrix language Mata
Here I give only a very brief introduction, for much more see:
William W. Gould (2018) The Mata Book: A Book for Serious Programmers
and Those Who Want to Be, College Station, TX: Stata Press.
You enter Mata by typing mata in Stata, and you leave it by typing end
In it you can work with matrices
. mata
------------------------------------------------- mata (type end to exit) -----
: foo = 3
: foo
3
:
: foo = 1, 2 \ 3, 4
: foo
1 2
+---------+
1 | 1 2 |
2 | 3 4 |
+---------+
:
: bar = J(2,2,1)
: bar
[symmetric]
1 2
+---------+
1 | 1 |
2 | 1 1 |
+---------+
: foo + bar
1 2
+---------+
1 | 2 3 |
2 | 4 5 |
+---------+
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- Mata
-------------------------------------------------------------------------------
Associative Array
However, matrices aren't the only thing that Mata can work with.
An ABM consist of agents that have properties, maybe a grid on which
agents live, etc.
A lot of the work in ABMs is just storing those properties, reading them,
and changing them.
We could do that in matrices, but associative arrays tend to be more
flexible
In an associative array you store key and value pairs. In the example
below the key is the last name, while the value is the first name. The
value could be anything.
. mata
------------------------------------------------- mata (type end to exit) -----
: name = AssociativeArray()
: name.put("Cox", "Nicholas")
: name.put("Jann", "Ben")
: name.put("Buis", "Maarten")
:
: name.get("Jann")
Ben
: end
-------------------------------------------------------------------------------
The key can also be a vector, leading to a 2, 3, or higher dimensional
array.
. mata
------------------------------------------------- mata (type end to exit) -----
: name.reinit("string",2)
: name.put(("Cox","middle"),"J.")
: name.put(("Cox", "first"), "Nicholas")
: name.put(("Jann", "first"), "Ben")
: name.put(("Buis", "middle"), "Leendert")
: name.put(("Buis", "first"), "Maarten")
:
: name.get(("Cox", "middle"))
J.
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- Mata
-------------------------------------------------------------------------------
Writing functions
When we write a new function we need to tell Mata what "kind of thing" it
can expect that function to return, and the kind of things the function
can expect as inputs
Things (scalars, matrices) we create inside a function are local to that
function
It is good practice to also first declare what kind of thing it is before
creating it.
. clear mata
.
. mata:
------------------------------------------------- mata (type end to exit) -----
:
: real scalar sum2(real scalar a, real scalar b)
> {
> real scalar result
>
> result = a + b
>
> return(result)
> }
:
: sum2(2,4)
6
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- Mata
-------------------------------------------------------------------------------
Loops, if
Looping and conditional execution are basic actions happing sometime in
any program
To illustrate the syntax used in Mata I add all even numbers between 1
and 10
. mata
------------------------------------------------- mata (type end to exit) -----
: result = 0
: for(i=1 ; i<=10; i++) {
> if (mod(i,2)==0) {
> result = result + i
> }
> }
: result
30
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- Mata
-------------------------------------------------------------------------------
structure
If you want to implement a larger project, like an ABM, then you will
want to split it up in sub-routines rather than put it all in one big
program.
However, you will often have quite a lot of "raw material" that you want
to share between sub-routines
In that case it is convenient to use a structure to organize that raw
material on a "tray", and pass that tray to the sub-routines.
. clear mata
.
. mata:
------------------------------------------------- mata (type end to exit) -----
: // ------------- define what goes where on our "tray"
: struct datastruct {
> real scalar a, b
> }
:
: // ------------------------------------ main function
: void arithmatic(real scalar a, real scalar b)
> {
> // create the tray called data
> struct datastruct scalar data
>
> // fill that tray
> data.a = a
> data.b = b
>
> // pass that tray on to sub-routines
> sum2(data)
> prod2(data)
> }
:
: // --------------------------------- the sub-routines
: real scalar sum2(struct datastruct scalar data)
> {
> return(data.a + data.b)
> }
:
: real scalar prod2(struct datastruct scalar data)
> {
> return(data.a * data.b)
> }
:
: // ---------------------------- use the main function
: arithmatic(2,4)
6
8
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- Mata
-------------------------------------------------------------------------------
Class
The next step is to use a class, which you can think of as a structure
that can also contain functions
The functions are local to the class
Each function has access to all material stored in the class
An "outsider", the user, can only interact with the public functions,
protected and private functions and variables can only be accessed from
within the class.
. clear mata
.
. mata:
------------------------------------------------- mata (type end to exit) -----
:
: class arithmatic {
> protected:
> real scalar a
> real scalar b
> real scalar sum()
> real scalar prod()
>
> public:
> void input()
> void run()
> }
:
: void arithmatic::input(real scalar aval, real scalar bval)
> {
> a = aval
> b = bval
> }
:
: real scalar arithmatic::sum()
> {
> return(a + b)
> }
:
: real scalar arithmatic::prod()
> {
> return(a * b)
> }
:
: void arithmatic::run()
> {
> sum()
> prod()
> }
:
: // how the user interacts with this program
: math = arithmatic()
: math.input(2,4)
: math.run()
6
8
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- Mata
-------------------------------------------------------------------------------
Inheritance
Classes can inherrit from other classes: All variables and functions in
the parent class also exist in the child class
The child class can add new functions
Function in the child class with the same name hide the function from the
parent.
That is a model for distributing common functions in ABMs
. mata
------------------------------------------------- mata (type end to exit) -----
: class arith extends arithmatic
> {
> private:
> real scalar ratio()
>
> public:
> void run()
> }
:
: real scalar arith::ratio()
> {
> return(a/b)
> }
:
: void arith::run() {
> sum()
> prod()
> ratio()
> }
:
: math = arith()
: math.input(2, 4)
: math.run()
6
8
.5
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- The abm_grid class
-------------------------------------------------------------------------------
Setting up a grid and finding neighbours
This class helps create and use a grid
The abm_grid class is available from GitHub:
https://github.com/maartenteaches/abm_grid
The main file is abm_grid.mata, all other files there are help files.
abm_grid
Below is an example that uses abm_grid to create a 10 by 10 grid that is
constant over time and allows for 1 agent per cell
. clear mata
. run abm_grid.mata
.
. mata
------------------------------------------------- mata (type end to exit) -----
: grid = abm_grid()
: grid.rdim(10)
: grid.cdim(10)
: grid.setup()
: end
-------------------------------------------------------------------------------
We use a grid in ABMs because it defines who the neighbours of an agent
are.
. mata
------------------------------------------------- mata (type end to exit) -----
: grid.find_ring(3,3,1)
1 2
+---------+
1 | 4 4 |
2 | 3 4 |
3 | 2 4 |
4 | 2 3 |
5 | 2 2 |
6 | 3 2 |
7 | 4 2 |
8 | 4 3 |
+---------+
: grid.find_ring(3,1,1)
1 2
+---------+
1 | 4 2 |
2 | 3 2 |
3 | 2 2 |
4 | 2 1 |
5 | 4 1 |
+---------+
: end
-------------------------------------------------------------------------------
| 1 2 3 4 | 1 2 3 4
--+------------ --+------------
1 | 1 |
2 | 5 4 3 2 | 4 3
3 | 6 o 2 3 | o 2
4 | 7 8 1 4 | 5 1
By default the cells that are outside the grid are ignored, but that can
be changed to let the grid "wrap around".
. mata
------------------------------------------------- mata (type end to exit) -----
: grid.torus(1)
: grid.setup()
:
: grid.find_ring(3,1,1)
1 2
+-----------+
1 | 4 2 |
2 | 3 2 |
3 | 2 2 |
4 | 2 1 |
5 | 2 10 |
6 | 3 10 |
7 | 4 10 |
8 | 4 1 |
+-----------+
: end
-------------------------------------------------------------------------------
| 1 2 3 4 .. 10
--+-----------------
1 |
2 | 4 3 5
3 | o 2 6
4 | 8 1 7
By default a neighbourhood with radius 1 is defined as all cells 1 step
removed, where a step could be horizontal, vertical, or diagonal. We
could change that to allow only horizontal and vertical steps.
. mata
------------------------------------------------- mata (type end to exit) -----
: grid.neumann(1)
: grid.setup()
:
: grid.find_ring(3,3,1)
1 2
+---------+
1 | 4 3 |
2 | 3 4 |
3 | 2 3 |
4 | 3 2 |
+---------+
: end
-------------------------------------------------------------------------------
| 1 2 3 4
--+------------
1 |
2 | 3
3 | 4 o 2
4 | 1
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- The abm_grid class
-------------------------------------------------------------------------------
Managing agents
In the grid we only store the agent id, the properties of those agents
(e.g. her color for Schelling seggregation model) should be stored in a
seperate associative array
The abm_grid class has function to check if an agent exists on a give
cell, to create an agent on a cell, and to retrieve the agent id of the
agent living at a cell.
. mata
------------------------------------------------- mata (type end to exit) -----
: grid.has_agent(1,1)
0
: grid.create_agent(1,1,5)
: grid.has_agent(1,1)
1
: grid.agent_id(1,1)
5
: end
-------------------------------------------------------------------------------
The abm_grid class also has functions to move, copy, or delete agents
. mata
------------------------------------------------- mata (type end to exit) -----
: //move agents
: grid.move_agent(1,1,2,2)
: grid.has_agent(1,1)
0
: grid.has_agent(2,2)
1
: grid.agent_id(2,2)
5
:
: //copy agents
: grid.copy_agent(2,2,1,1)
: grid.agent_id(2,2)
5
: grid.agent_id(1,1)
5
:
: // delete agents
: grid.delete_agent(2,2)
: grid.has_agent(2,2)
0
: end
-------------------------------------------------------------------------------
After the ABM has run we often want to export the results to Stata to
create graphs
. mata
------------------------------------------------- mata (type end to exit) -----
: grid.create_agent(5,3,2)
: grid.create_agent(2,9,1)
: grid.extract()
1 2 3 4 5
+---------------------+
1 | 1 1 0 1 5 |
2 | 2 9 0 1 1 |
3 | 5 3 0 1 2 |
+---------------------+
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- The abm_grid class
-------------------------------------------------------------------------------
Adding a time dimension to the grid
We need to add a time dimension to our grid if agents can move around on
the grid, e.g. in Schelling's segregation model.
. mata
------------------------------------------------- mata (type end to exit) -----
: grid.tdim(5)
: grid.setup()
:
: // we setup the intial grid
: grid.create_agent(1,1,5, 0)
: grid.create_agent(3,4,1, 0)
: grid.create_agent(4,4,2, 0)
:
: // next iteration we first copy the grid
: grid.copy_grid(0,1)
: // and make changes as necessary
: grid.move_agent(1,1,2,3,1)
:
: //at the end we can recover the entire history of the grid:
: grid.extract()
1 2 3 4 5
+---------------------+
1 | 1 1 0 1 5 |
2 | 2 3 1 1 5 |
3 | 3 4 0 1 1 |
4 | 3 4 1 1 1 |
5 | 4 4 0 1 2 |
6 | 4 4 1 1 2 |
+---------------------+
: end
-------------------------------------------------------------------------------
If the agents cannot move, e.g. in the immunization model, then grid does
not need a time dimension.
Now the associative array storing the agents needs to have a time
dimension.
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- The abm_grid class
-------------------------------------------------------------------------------
Multiple agents per cell
For some ABMs we want multiple agents to live in the same cell
For example, a predator-prey model where sheep and wolve agents wander on
a grid
Each agent still needs a distinct "adress", so we can add a "location"
dimension within cells
. mata
------------------------------------------------- mata (type end to exit) -----
: grid.idim(3)
: grid.setup()
:
: // add three agents to the cell 1,1 at itereation 0
: grid.create_agent(1,1,5, 0,1)
: grid.create_agent(1,1,6, 0,2)
: grid.create_agent(1,1,7, 0,3)
:
: grid.agent_id(1,1,0,1)
5
: grid.agent_ids(1,1,0)
1 2 3
+-------------+
1 | 5 6 7 |
+-------------+
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Introduction -- The abm_grid class
-------------------------------------------------------------------------------
How do I use it?
A typical start of an ABM would be:
. clear mata
. run abm_grid.mata
.
. mata
------------------------------------------------- mata (type end to exit) -----
: class grid extends abm_grid
> {
>
> }
:
: class abm
> {
> private:
> class grid scalar universe
>
> // other things required for the model
> }
: end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
How to start
I start with a rough idea of what the model should be, e.g.:
In this model there are blue and red agents, that live in cells on a
square grid. Each cell can contain only one agent, but not all cells
are occupied. An agent is happy if the proportion of agents of the
same colour in her neighbourhood is above a certain cut-off value. If
she is unhappy she will move to the nearest available cell that does
satisfy her needs.
From that I start defining the class, what kind of variables do I need to
store, what kind of functions do I need.
This class definition becomes the todo list for the project.
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
void rdim()
void cdim()
void tdim()
void limit()
void nblue()
void nred()
// run the model
void run()
// some functions to export or summarize the results
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
Receiving and returning settings
We want the user to be able to specify
the number of rows and columns in the grid
the number of iterations
the number of red and blue agents
the limit, below which the agent becomes unhappy
Lets start with the number of rows:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
void rdim()
void cdim()
void tdim()
void limit()
void nblue()
void nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------------ setup grid
void schelling::rdim(real scalar dim)
{
string scalar errmsg
if (dim <= 0 | mod(dim,1)!= 0 ){
errmsg = "dim must be a positive integer"
_error(3300, errmsg)
}
universe.rdim(dim)
rdim = dim
}
end
We are checking the input to see if it is a positive integer, and we will
have to do that for most of the other settings as well. Lets put that in
a separate function:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
void rdim()
void cdim()
void tdim()
void limit()
void nblue()
void nred()
// run the model
void run()
// some functions to export or summarize the results
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
void schelling::rdim(real scalar dim)
{
is_posint(dim, "dim")
universe.rdim(dim)
rdim = dim
}
end
Can't we make it such that if we give rdim() an argument it sets rdim to
that value, but we don't give rdim() an argument it returns the value of
rdim?
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
schel = schelling()
schel.rdim(5)
schel.rdim()
end
That looks good, just extend to all other settings
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
Setting up
Seting up the universe is something we can let abm_grid do, the only
thing we need to add is to make sure tdim is specified:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
end
To setup the agents we first need to check some of the options:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
}
end
Now we need to change the agents associative array such that the key is
the agent id (a real number) and fill it with red and blue agents:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
end
The final step is put these agents on random places on the grid in the
initial iteration:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
Summarize the neighbourhood
Now we agents living on the grid. The next step would be for them to look
around and see if they are happy with their neighbourhood.
Lets start with a routine that returns the color of an agent living at a
certain cell:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
Next we loop over all neighbours and count the proportion of same color
neighbours:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c)
{
real scalar k, s, i, j, l
real matrix neigh
string scalar self
self = get_color(t,r,c)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
if (universe.has_agent(i,j,t) ) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return(s/k)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
Right now a missing value will be returned if an agent has no neighbours,
we need to fix that :
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c)
{
real scalar k, s, i, j, l
real matrix neigh
string scalar self
self = get_color(t,r,c)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
if (universe.has_agent(i,j,t) ) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
We are later also going to use this function to find new locations. In
that case the focal cell is not the cell where the agent currently lives:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
if (universe.has_agent(i,j,t) ) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
In that case the neighbourhood of the focal cell could contain the cell
where the agent currently lives. We would want to exclude that cell, as
she will be moving:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
Now we can create another function that runs sum_neighbourhood() and
compares the returned proportion with the limit:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
Move to a better cell
There is an function in abm_grid that can move an agent from one location
to another, so we only need to focus on finding the nearest location that
will make the agent happy.
Lets start with a loop that moves over all cells in the grid in an
outward spiral:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius
real matrix neigh
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
}
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
Now we can check any empty cells to see if the agent will be happy there,
and move when that is the case:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd
real matrix neigh
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
}
}
}
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
We can stop the loop the first time we find a appropriate cell:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
Run the simulation
First we will define the action of a single agent at each itteration:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
Next we will define the function run() that lets the user start the
simulation. It starts by setting up the universe and agents, populate
that universe:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
setup_universe()
setup_agents()
populate_universe()
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
Than it loops over iterations, each time copying the grid to the next
iteration, and than looping over the cells of that grid:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
If there is an agent in that cell it lets it take its step:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
// some functions to export or summarize the results
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
display summary statistics
We want to display a table that shows for each iteration the average
proportion of same color neighbours and the proportion of happy agents.
Start with creating a matrix with the first column indicating the
iteration:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prob()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
}
// -------------------------------------------------------------- return results
real matrix schelling::summ()
{
real matrix res
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0 // number of neighbours
s = 0 // number of same color neighbours
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
end
Next we fill in the second column with the average proportion of same
color agents:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
}
// -------------------------------------------------------------- return results
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,1,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k] = sum_neighbours(t,r,c,t,r,c)
k = k+1
}
}
}
res[t+1, 2] = mean(prop)
}
return(res)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summ()
end
Now we can fill the last column:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
}
// -------------------------------------------------------------- return results
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summ()
end
If export this matrix to Stata then we can use matlist to make this table
look nice:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
void summtable()
}
// -------------------------------------------------------------- return results
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summtable()
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
Export to Stata
To plot the agents on the grid we first need to export it to Stata. We
would like for each cell with an agent the:
row number
column number
iteration number
id of that agent
the color of that agent
whether or not the agent is happy
So this function must receive 6 names of variables that are to be created
With the abm_grid class we can extract the row, column, iteration number
and the id, though the 4{sup:th} column is superfluous: That contains the
location within the cell, but we only allowed one agent per cell
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
void summtable()
void extract()
}
// -------------------------------------------------------------- return results
void schelling::extract(string rowvector names)
{
real matrix res
if(cols(names)!=6){
_error(3300, "names must contain 6 columns")
}
res = universe.extract()
res = res[.,1..3], res[.,5]
}
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summtable()
end
The 5{sup:th} column is the id of the agent, which we can use to recover
the color from the associative array agents:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
void summtable()
void extract()
}
// -------------------------------------------------------------- return results
void schelling::extract(string rowvector names)
{
real matrix res
real colvector color
real scalar id, i
if(cols(names)!=6){
_error(3300, "names must contain 6 columns")
}
res = universe.extract()
color = J(rows(res),1,.)
for (i= 1 ; i <= rows(res) ; i++) {
id = res[i,5]
color[i] = ( agents.get(id) == "red" )
}
res = res[.,1..3], res[.,5], color
}
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summtable()
end
We can now add the whether or not the agent is happy:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
void summtable()
void extract()
}
// -------------------------------------------------------------- return results
void schelling::extract(string rowvector names)
{
real matrix res
real colvector color, happy
real scalar id, i, t, r, c
if(cols(names)!=6){
_error(3300, "names must contain 6 columns")
}
res = universe.extract()
color = J(rows(res),1,.)
happy = J(rows(res),1,.)
for (i= 1 ; i <= rows(res) ; i++) {
id = res[i,5]
color[i] = ( agents.get(id) == "red" )
r = res[i,1]
c = res[i,2]
t = res[i,3]
happy[i] = is_happy(t,r,c,t,r,c)
}
res = res[.,1..3], res[.,5], color, happy
}
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summtable()
end
Now we need to store that matrix as a Stata dataset:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
void summtable()
void extract()
}
// -------------------------------------------------------------- return results
void schelling::extract(string rowvector names)
{
real matrix res
real colvector color, happy
real rowvector idx
real scalar id, i, t, r, c
if(cols(names)!=6){
_error(3300, "names must contain 6 columns")
}
res = universe.extract()
color = J(rows(res),1,.)
happy = J(rows(res),1,.)
for (i= 1 ; i <= rows(res) ; i++) {
id = res[i,5]
color[i] = ( agents.get(id) == "red" )
r = res[i,1]
c = res[i,2]
t = res[i,3]
happy[i] = is_happy(t,r,c,t,r,c)
}
res = res[.,1..3], res[.,5], color, happy
idx = st_addvar("float", names)
if (st_nobs() < rows(res) ) {
st_addobs(rows(res) - st_nobs())
}
st_store(.,idx,res)
}
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summtable()
end
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- the Schellings segregation model
-------------------------------------------------------------------------------
Plot the grid
We can use {cmd twoway scatter} to plot the grid for the initial
iteration:
clear mata
drop _all
run abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours()
real scalar is_happy()
void move()
void step()
public:
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
void run()
void extract()
real matrix summ()
void summtable()
}
// -------------------------------------------------------------- error messages
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar k, i
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
// --------------------------------------------------------------- run the model
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, stop, rd, cd, maxradius
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1,) == 0) {
if( is_happy(t+1,rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
// ---------------------------------------------------------------------- export
void schelling::extract(string rowvector names)
{
real matrix res
real colvector color, happy
real rowvector idx
real scalar i, k, t,r,c
if(cols(names)!=6){
_error(3300, "names must contain 6 columns")
}
res = universe.extract()
color = J(rows(res),1,.)
happy = J(rows(res),1,.)
k = cols(res)
for (i = 1 ; i <= rows(res); i++) {
color[i] = agents.get(res[i,k]) == "red"
t = res[i,3]
r = res[i,1]
c = res[i,2]
happy[i] = is_happy(t,r,c,t,r,c)
}
res = res[.,1..3], res[.,5], color, happy
idx = st_addvar("float", names)
if (st_nobs() < rows(res) ) {
st_addobs(rows(res) - st_nobs())
}
st_store(.,idx,res)
}
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, r, c, k
res = J(tdim+1,3,.)
res[.,1] = (0 .. tdim)'
for(t=0 ; t<= tdim; t++) {
prop = J(nagents, 2, .)
k = 1
for(r=1; r <= rdim; r++) {
for (c=1 ; c<=cdim ; c++ ) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c, t, r, c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k + 1
}
}
}
res[|t+1,2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
end
There are still some improvements necessary:
We don't need to label the axes
It is a square grid, so the graph should also be square
We could add lines to highlight the cells
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
void summtable()
void extract()
}
// -------------------------------------------------------------- return results
void schelling::extract(string rowvector names)
{
real matrix res
real colvector color, happy
real rowvector idx
real scalar id, i, t, r, c
if(cols(names)!=6){
_error(3300, "names must contain 6 columns")
}
res = universe.extract()
color = J(rows(res),1,.)
happy = J(rows(res),1,.)
for (i= 1 ; i <= rows(res) ; i++) {
id = res[i,5]
color[i] = ( agents.get(id) == "red" )
r = res[i,1]
c = res[i,2]
t = res[i,3]
happy[i] = is_happy(t,r,c,t,r,c)
}
res = res[.,1..3], res[.,5], color, happy
idx = st_addvar("float", names)
if (st_nobs() < rows(res) ) {
st_addobs(rows(res) - st_nobs())
}
st_store(.,idx,res)
}
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summtable()
foo.extract(("r","c","t", "id", "red", "happy"))
end
replace r = -r
twoway scatter r c if red == 0 & t == 0, ///
msymbol(O) mcolor(blue) || ///
scatter r c if red == 1 & t == 0, ///
msymbol(O) mcolor(red) ///
yscale(off range(-10.5 -.5)) ///
xscale(off range(.5 10.5)) ///
ylab(none) xlab(none) ///
plotregion(margin(zero)) ///
aspect(1) ///
xline(.5(1)9.5) ///
yline(-.5(-1)-9.5) ///
legend(off)
We can improve the graph further:
We can distinguish between happy and non-happy agents
We can add the agent id
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
void summtable()
void extract()
}
// -------------------------------------------------------------- return results
void schelling::extract(string rowvector names)
{
real matrix res
real colvector color, happy
real rowvector idx
real scalar id, i, t, r, c
if(cols(names)!=6){
_error(3300, "names must contain 6 columns")
}
res = universe.extract()
color = J(rows(res),1,.)
happy = J(rows(res),1,.)
for (i= 1 ; i <= rows(res) ; i++) {
id = res[i,5]
color[i] = ( agents.get(id) == "red" )
r = res[i,1]
c = res[i,2]
t = res[i,3]
happy[i] = is_happy(t,r,c,t,r,c)
}
res = res[.,1..3], res[.,5], color, happy
idx = st_addvar("float", names)
if (st_nobs() < rows(res) ) {
st_addobs(rows(res) - st_nobs())
}
st_store(.,idx,res)
}
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summtable()
foo.extract(("r","c","t", "id", "red", "happy"))
end
replace r = -r
twoway scatter r c if red == 0 & t == 0 & happy == 1, ///
msymbol(O) mcolor(blue) msize(*2) ///
mlabel(id) mlabpos(0) || ///
scatter r c if red == 1 & t == 0 & happy == 1, ///
msymbol(O) mcolor(red) msize(*2) ///
mlabel(id) mlabpos(0) || ///
scatter r c if red == 0 & t == 0 & happy == 0, ///
msymbol(Oh) mcolor(blue) msize(*2) ///
mlabel(id) mlabpos(0) || ///
scatter r c if red == 1 & t == 0 & happy == 0, ///
msymbol(Oh) mcolor(red) msize(*2) ///
mlabel(id) mlabpos(0) ///
yscale(off range(-10.5 -.5)) ///
xscale(off range(.5 10.5)) ///
ylab(none) xlab(none) ///
plotregion(margin(zero)) ///
aspect(1) ///
xline(.5(1)9.5) ///
yline(-.5(-1)-9.5) ///
legend(off)
Now we can use a loop to create the same graph for all iterations:
clear all
do abm_grid.mata
mata
mata set matastrict on
class grid extends abm_grid
{
}
class schelling {
private:
// all variables are private
class grid scalar universe
class AssociativeArray scalar agents
real scalar rdim
real scalar cdim
real scalar tdim
real scalar limit
real scalar nblue
real scalar nred
real scalar nagents
void is_posint()
void is_prop()
// initialize the model
void setup_agents()
void setup_universe()
void populate_universe()
string scalar get_color()
real scalar sum_neighbours() // proportion of same color in neighbourhood
real scalar is_happy() // is that proportion larger than limit
void move() // search for nearest empty spot where agent would be happy
void step() // a single step by a single agent
public:
// set or return settings
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic limit()
transmorphic nblue()
transmorphic nred()
// run the model
void run()
real matrix summ()
void summtable()
void extract()
}
// -------------------------------------------------------------- return results
void schelling::extract(string rowvector names)
{
real matrix res
real colvector color, happy
real rowvector idx
real scalar id, i, t, r, c
if(cols(names)!=6){
_error(3300, "names must contain 6 columns")
}
res = universe.extract()
color = J(rows(res),1,.)
happy = J(rows(res),1,.)
for (i= 1 ; i <= rows(res) ; i++) {
id = res[i,5]
color[i] = ( agents.get(id) == "red" )
r = res[i,1]
c = res[i,2]
t = res[i,3]
happy[i] = is_happy(t,r,c,t,r,c)
}
res = res[.,1..3], res[.,5], color, happy
idx = st_addvar("float", names)
if (st_nobs() < rows(res) ) {
st_addobs(rows(res) - st_nobs())
}
st_store(.,idx,res)
}
real matrix schelling::summ()
{
real matrix res, prop
real scalar t, k, r, c
res = J(tdim+1, 3, .)
res[.,1] = (0..tdim)'
for(t=0; t <= tdim ; t++){
prop = J(nagents,2,.)
k=1
for (r=1; r<=rdim ; r++) {
for (c=1 ; c <= cdim ; c++) {
if (universe.has_agent(r,c,t)) {
prop[k,1] = sum_neighbours(t,r,c,t,r,c)
prop[k,2] = is_happy(t,r,c,t,r,c)
k = k+1
}
}
}
res[|t+1, 2 \ t+1,3|] = mean(prop)
}
return(res)
}
void schelling::summtable(| string scalar matname)
{
real matrix res
string scalar cmd
res = summ()
if (args()==0) {
matname = st_tempname()
}
st_matrix(matname,res)
st_matrixcolstripe(matname, (J(3,1, "") , ("Iteration" \ "prop_same" \ "prop_happy")) )
cmd = "matlist " + matname + ", underscore format(%10.3g) names(columns) " +
"border(rows) " +
"title(Average proportion of neighbours with same color and average proportion of happy agents)"
stata(cmd)
}
// --------------------------------------------------------------- run the model
void schelling::step(real scalar t, real scalar r, real scalar c)
{
if (is_happy(t,r,c,t,r,c)==0) {
move(t,r,c)
}
}
void schelling::run() {
real scalar t, k, r, c
real matrix coords
setup_universe()
setup_agents()
populate_universe()
for(t=0; t < universe.tdim() ; t++) {
universe.copy_grid(t, t+1)
coords = universe.schedule()
for(k=1; k <= rows(coords); k++){
r = coords[k,1]
c = coords[k,2]
if (universe.has_agent(r,c,t)) {
step(t, r, c)
}
}
}
}
// ------------------------------------------------------------ helper functions
void schelling::move(real scalar t, real scalar r, real scalar c)
{
real scalar radius, i, maxradius, rd, cd, stop
real matrix neigh
stop = 0
maxradius = max( ( rdim-r, r-1, cdim-c, c-1 ) )
for (radius = 1 ; radius <= maxradius ; radius ++) {
neigh = universe.find_ring(r,c,radius)
for (i = 1; i <= rows(neigh); i++) {
rd = neigh[i,1]
cd = neigh[i,2]
if (universe.has_agent(rd,cd, t+1) == 0) {
if ( is_happy(t+1, rd, cd, t, r, c) ) {
universe.move_agent(r,c, rd, cd, t+1)
stop = 1
break
}
}
}
if (stop == 1) break
}
}
string scalar schelling::get_color(real scalar t, real scalar r, real scalar c)
{
real scalar id
if (universe.has_agent(r,c,t)) {
id = universe.agent_id(r, c,t)
}
else {
_error(3300, "no agent in this cell")
}
return(agents.get(id))
}
real scalar schelling::sum_neighbours(
real scalar t, real scalar r, real scalar c,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar k, s, i, j, l, isself
real matrix neigh
string scalar self
self = get_color(ts,rs,cs)
neigh = universe.find_ring(r,c,1)
k = 0
s = 0
for (l=1 ; l <=rows(neigh) ; l++) {
i = neigh[l,1]
j = neigh[l,2]
isself = (i==rs & j==cs)
if (universe.has_agent(i,j,t) & !isself) {
k = k + 1
s = s + (get_color(t,i,j)==self)
}
}
return((k==0 ? 1 : s/k))
}
real scalar schelling::is_happy(
real scalar t , real scalar r , real scalar c ,
real scalar ts, real scalar rs, real scalar cs)
{
real scalar prop
prop = sum_neighbours(t,r,c,ts,rs,cs)
return(prop >= limit)
}
// -------------------------------------------------------------- error checking
void schelling::is_posint(real scalar val, string scalar name)
{
string scalar errmsg
if (val <= 0 | mod(val,1)!= 0 ){
errmsg = name + " must be a positive integer"
_error(3300, errmsg)
}
}
void schelling::is_prop(real scalar val, string scalar name)
{
string scalar errmsg
if (val < 0 | val > 1 ){
errmsg = name + " must be between or including 0 and 1"
_error(3300, errmsg)
}
}
// ------------------------------------------------------------------ setup grid
transmorphic schelling::rdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.rdim(dim)
rdim = dim
}
else {
return(universe.rdim())
}
}
transmorphic schelling::cdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.cdim(dim)
cdim = dim
}
else {
return(universe.cdim())
}
}
transmorphic schelling::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim,"dim")
universe.tdim(dim)
tdim = dim
}
else {
return(universe.tdim())
}
}
void schelling::setup_universe()
{
if (tdim()==.) {
_error(3000, "tdim must be set first")
}
universe.setup()
}
// ---------------------------------------------------------------- setup agents
transmorphic schelling::limit(| real scalar val)
{
if (args() == 1) {
is_prop(val,"value")
limit = val
}
else {
return(limit)
}
}
transmorphic schelling::nblue(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nblue = val
}
else {
return(nblue)
}
}
transmorphic schelling::nred(| real scalar val)
{
if (args() == 1) {
is_posint(val,"value")
nred = val
}
else {
return(nred)
}
}
void schelling::setup_agents()
{
real scalar i, k
string scalar errmsg
if (limit == .) {
_error(3300, "limit has not been set")
}
if (nblue == .) {
_error(3300, "number of blue agents has not been set")
}
if (nred == .) {
_error(3300, "number of red agents has not been set")
}
nagents = nred + nblue
if (nagents > universe.size()) {
errmsg = "you specified " + strofreal(nagents) + " agents on a " +
strofreal(rdim()) + " by " + strofreal(cdim()) + " grid"
_error(3300, errmsg)
}
agents.reinit("real")
k = 0
for (i = 1 ; i <= nred ; i++) {
k = k +1
agents.put(k, "red")
}
for (i = 1 ; i <= nblue ; i++) {
k = k +1
agents.put(k, "blue")
}
}
void schelling::populate_universe(){
real matrix grid0
real scalar i, r, c
grid0 = universe.schedule()
_jumble(grid0)
for(i = 1 ; i <= nagents ; i++) {
r = grid0[i,1]
c = grid0[i,2]
universe.create_agent(r,c,i,0)
}
}
foo = schelling()
foo.rdim(10)
foo.cdim(10)
foo.tdim(5)
foo.nred(45)
foo.nblue(45)
foo.limit(.3)
foo.run()
foo.summtable()
foo.extract(("r","c","t", "id", "red", "happy"))
end
replace r = -r
forvalues t = 0/5 {
twoway scatter r c if red == 0 & t == `t' & happy == 1, ///
msymbol(O) mcolor(blue) msize(*2) ///
mlabel(id) mlabpos(0) || ///
scatter r c if red == 1 & t == `t' & happy == 1, ///
msymbol(O) mcolor(red) msize(*2) ///
mlabel(id) mlabpos(0) || ///
scatter r c if red == 0 & t == `t' & happy == 0, ///
msymbol(Oh) mcolor(blue) msize(*2) ///
mlabel(id) mlabpos(0) || ///
scatter r c if red == 1 & t == `t' & happy == 0, ///
msymbol(Oh) mcolor(red) msize(*2) ///
mlabel(id) mlabpos(0) ///
yscale(off range(-10.5 -.5)) ///
xscale(off range(.5 10.5)) ///
ylab(none) xlab(none) ///
plotregion(margin(zero)) ///
aspect(1) ///
xline(.5(1)9.5) ///
yline(-.5(-1)-9.5) ///
legend(off) ///
name(gr`t', replace) title(Iteration `t')
}
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Building Agent Based Models -- Other models
-------------------------------------------------------------------------------
Other models
The code for the immunization model is here:
set seed 123456789
set rmsg on
clear mata
drop _all
// -----------------------------------class definitions
// import abm_grid class definition
do abm_grid.mata
local vulnerable = 0
local infected = 1
local infectious = 2
local healed = 3
// this allows you to change or add to the grid class when needed
mata
class grid extends abm_grid
{
}
// the main agent based model
class abm
{
private:
class grid scalar universe
class AssociativeArray scalar agents
real scalar not_vaccinated
real scalar tdim
real scalar radius
real matrix infected
real scalar n_infected
void setup_grid()
void setup_agents()
void create_agents()
void populate_universe()
void copy_agents()
void infect()
real scalar status()
real scalar is_vaccinated()
void initial_infection()
public:
transmorphic rdim()
transmorphic cdim()
transmorphic tdim()
transmorphic torus()
transmorphic not_vaccinated()
transmorphic radius()
void run()
real scalar n_infected()
void extract()
}
// --------------------------------------setup the grid
transmorphic abm::rdim(| real scalar dim)
{
if (args() == 1) {
universe.rdim(dim)
}
else {
return(universe.rdim())
}
}
transmorphic abm::cdim(| real scalar dim)
{
if (args() == 1) {
universe.cdim(dim)
}
else {
return(universe.cdim())
}
}
transmorphic abm::torus(| real scalar val)
{
if (args() == 1) {
universe.torus(val)
}
else {
return(universe.torus())
}
}
void abm::setup_grid()
{
if (universe.torus() == . ) {
universe.torus(1)
}
universe.setup()
}
// ------------------------------------- setup agents
transmorphic abm::not_vaccinated(| real scalar dim)
{
if (args() == 1) {
if (dim <= 0 | mod(dim,1) != 0 ) {
_error(3000, "dim must be an integer larger than 0")
}
not_vaccinated = dim
}
else {
return(not_vaccinated)
}
}
transmorphic abm::tdim(| real scalar dim)
{
if (args() == 1) {
if (dim <= 0 | mod(dim,1) != 0) {
_error(3000, "dim must be an integer larger than 0")
}
tdim = dim
}
else {
return(tdim)
}
}
transmorphic abm::radius(| real scalar number)
{
if (args() == 1) {
if (number <= 0 | mod(number,1) != 0) {
_error(3000, "number must be an integer larger than 0")
}
radius = number
}
else {
return(radius)
}
}
void abm::setup_agents()
{
if (tdim == .) {
_error(3000, "tdim must be specified")
}
if (not_vaccinated == .) {
_error(3000, "not_vaccinated must be specified")
}
if (not_vaccinated > universe.size()) {
_error(3000, "specified more not vaccinated agents than cells on grid")
}
if (radius == .) {
radius = 2
}
agents.reinit("real", 2)
}
// -------------------------------- helper functions
real scalar abm::status(real scalar r, real scalar c, real scalar t)
{
real scalar id
real rowvector key
if (universe.has_agent(r,c)) {
id = universe.agent_id(r,c)
key = id, t
return(agents.get(key))
}
else {
return(0)
}
}
real scalar abm::is_vaccinated(real scalar r, real scalar c, real scalar t)
{
return(!universe.has_agent(r,c))
}
void abm::copy_agents(real scalar t)
{
real scalar i, content
real rowvector key
for(i=1 ; i <= not_vaccinated ; i++) {
key = i,t
content = agents.get(key)
key = i,t+1
agents.put(key,content)
}
}
real scalar abm::n_infected() {
return(n_infected)
}
void abm::extract(string rowvector names)
{
real matrix grid, ag
real rowvector idx
real colvector status
real scalar i
if (cols(names) != 5) {
_error(3000, "5 names need to be specified")
}
grid = universe.extract()
_sort(grid,5)
grid = grid[.,(1,2)]
grid = J(tdim+1,1,grid)
ag = agents.keys()
_sort(ag,(2,1))
status = J(rows(ag),1,.)
for(i=1 ; i<= rows(ag); i++) {
status[i] = agents.get(ag[i,.])
}
ag = ag, status, grid
idx = st_addvar("float", names)
if (st_nobs() < rows(ag) ) {
st_addobs(rows(ag) - st_nobs())
}
st_store(.,idx,ag)
}
// ----------------------------------- run the model
void abm::create_agents()
{
real scalar i, id
real matrix unvac
for (i = 1 ; i <= not_vaccinated ; i++) {
agents.put((i,1),0)
}
}
void abm::populate_universe()
{
real scalar i
real matrix grid
grid = jumble(universe.schedule())
infected = grid[1,.]
for(i = 1 ; i <= not_vaccinated ; i++) {
universe.create_agent(grid[i,1],grid[i,2],i)
}
}
void abm::initial_infection()
{
real scalar id
id = universe.agent_id(infected[1], infected[2])
agents.put((id,1),`infected')
n_infected = 1
}
void abm::infect(real scalar t)
{
real scalar i, id, r, c, j, rn, cn, idn
real rowvector key
real colvector tokeep
real matrix neighbours, toadd
tokeep = J(rows(infected),1,1)
toadd = J(0,2,.)
for(i = 1 ; i <= rows(infected) ; i++) {
r = infected[i,1]
c = infected[i,2]
id = universe.agent_id(r, c)
if ( status(r, c, t) == `infectious' ) {
agents.put((id,t+1),`healed')
tokeep[i] = 0
}
else{
key = id, t+1
agents.put(key, `infectious')
neighbours = universe.find_spiral(r,c,radius)
for (j = 1 ; j <= rows(neighbours); j++) {
rn = neighbours[j,1]
cn = neighbours[j,2]
if ( is_vaccinated(rn,cn,t+1) == 0 & status(rn,cn,t+1) == 0 ) {
idn = universe.agent_id(rn,cn)
key = idn,t + 1
agents.put(key,`infected')
n_infected = n_infected + 1
toadd = toadd \ (rn,cn)
}
}
}
}
infected = select(infected, tokeep)
infected = infected \ toadd
}
void abm::run()
{
real scalar t
setup_grid()
setup_agents()
create_agents()
populate_universe()
initial_infection()
for(t=1; t <= tdim ; t++) {
copy_agents(t)
infect(t)
}
}
end
// run a single model and graph the evolution of the infection over time
mata
foo = abm()
foo.rdim(40)
foo.cdim(40)
foo.tdim(24)
foo.radius(2)
foo.not_vaccinated(240)
foo.run()
foo.extract(("id", "time", "status", "y", "x"))
end
replace y = - y
mata st_local("rdim", strofreal(foo.rdim()))
mata st_local("cdim", strofreal(foo.cdim()))
mata st_local("tdim", strofreal(foo.tdim()))
tempfile data
save `data'
drop _all
set obs `rdim'
gen x = _n
expand `cdim'
bys x : gen y = -_n
expand `tdim' + 1
bys y x : gen time = _n
merge 1:1 x y time using `data'
assert _merge != 2
drop _merge
forvalues t = 1/25 {
scatter y x if status == 0 & time == `t' , ///
msymbol(Oh) || ///
scatter y x if status == 1 & time == `t' , ///
msymbol(S) mcolor(gs8) || ///
scatter y x if status == 2 & time == `t' , ///
msymbol(S) || ///
scatter y x if status == 3 & time == `t' , ///
msymbol(+) || ///
scatter y x if status == . & time == `t', ///
msymbol(o) msize(*.10) mcolor(gs8) ///
aspect(1) name(gr`t', replace) ///
yscale(off range(-41.5 -0.5)) ///
xscale(off range(.5 41.5)) ///
plotregion(margin(zero)) ///
title("Iteration `t'") ///
legend(order( 1 "vulnerable" ///
2 "infected" ///
3 "infectious" ///
4 "healed/immunized" ///
5 "immunized"))
}
// run multiple models and plot the seriousness of the infection against
// the immunization rate
mata
scen = (16, 80, 120, 160, 200, 240, 280, 320, 360, 400)
res = J(200,2,.)
k=0
for (i = 1; i <= 10 ; i++) {
printf("number of infected: %3.0f \n", scen[i])
foo.not_vaccinated(scen[i])
for(j = 1; j<= 20; j++) {
k = k +1
foo.run()
res[k,.] = scen[i], foo.n_infected()
}
}
idx = st_addvar("float", ("N","n"))
st_addobs(200)
st_store((1..200)',idx,res)
end
gen perc_notvac = N/(40*40)*100
gen perc_inf = n/N*100
scatter perc_inf perc_notvac, msymbol(o) mcolor(%20) ///
xtitle("percentage not vaccinated") ///
ytitle("percentage infected of not vacinated ") ///
name(scenarios, replace)
The code for the bot and public mood model is here:
clear all
set rmsg on
local content = 1
local clicks = 2
mata
mata set matastrict on
class clicks {
private:
class AssociativeArray scalar agents
class AssociativeArray scalar items
class AssociativeArray scalar pr
real scalar Nagents
real scalar newscycle
real scalar pr_good
real scalar tdim
real scalar baseclicks
real scalar botclicks
void is_pr()
void is_posint()
void setup()
void populate_agents()
void populate_news()
void news_step()
void agent_step()
void pr()
public:
transmorphic Nagents()
transmorphic newscycle()
transmorphic pr_good()
transmorphic tdim()
transmorphic baseclicks()
transmorphic botclicks()
void run()
transmorphic mood()
transmorphic items()
}
transmorphic clicks::items(| string rowvector names )
{
real matrix res
real rowvector idx
real scalar i
res = items.keys()
res = res,J(rows(res), 1, .)
for (i=1 ; i <= rows(res); i++) {
res[i,4] = items.get(res[|i,1 \ i,3|])
}
_sort(res, (1,2,3))
if (args() == 1) {
if (cols(names) != 4) {
_error(3300, "4 names must be specified")
}
idx = st_addvar("float", names)
if (st_nobs() < rows(res) ) {
st_addobs(rows(res) - st_nobs())
}
st_store(.,idx,res)
}
else {
return(res)
}
}
transmorphic clicks::mood(| string rowvector names)
{
real scalar i, t, c
real rowvector key, idx
real matrix res
res = J(tdim+1,10,0)
for (t=0 ; t <= tdim; t++) {
for (i = 1 ; i<= Nagents; i++) {
key = i,t
c = agents.get(key)
res[t+1,c] = res[t+1,c] + 1
}
}
if (args()==1) {
if (cols(names) != 10) {
_error(3300, "10 names must be specified")
}
idx = st_addvar("float", names)
if (st_nobs() < rows(res) ) {
st_addobs(rows(res) - st_nobs())
}
st_store(.,idx,res)
}
else {
return(res)
}
}
void clicks::agent_step(real scalar i, real scalar t)
{
real scalar u, j, pref, prob
real rowvector key
key = i, t-1
pref = agents.get(key)
key = i, t
agents.put(key,pref)
u = runiform(1,1)
prob = 0
for(j=t+1; j <= t + newscycle; j++) {
prob = prob + pr.get(j)
if (u < prob) {
key = j, t, `content'
if (items.get(key)) {
key = i,t
agents.put(key,min((agents.get(key) +1, 10)))
}
else {
key = i,t
agents.put(key,max((agents.get(key)-1, 1)))
}
key = j, t, `clicks'
items.put(key,items.get(key)+1)
break
}
}
}
void clicks::news_step(real scalar t)
{
real scalar i, content, good
real rowvector key
for(i = t+1 ; i <= t + newscycle - 1; i++) {
key = i, t-1, `content'
content = items.get(key)
key = i, t, `content'
items.put(key, content)
key = i, t-1, `clicks'
content = items.get(key)
key = i, t, `clicks'
items.put(key,content)
}
key = t + newscycle, t, `content'
good = runiform(1,1) < pr_good
items.put(key, good)
key = t + newscycle, t, `clicks'
items.put(key, ceil(runiform(1,1)*baseclicks) + (1-good)*botclicks)
}
void clicks::pr(real scalar t)
{
real scalar i, n, n_clicks
real colvector key
n_clicks = 0
for(i = t+1 ; i <= t + newscycle ; i++) {
key = i,t,`clicks'
n = items.get(key)
n_clicks = n_clicks + n
}
for(i = t+1 ; i <= t + newscycle ; i++) {
key = i,t,`clicks'
n = items.get(key)
pr.put(i, n/n_clicks)
}
}
void clicks::is_posint(real scalar val, string scalar name, | string scalar zero_ok)
{
string scalar errmsg
if (args() == 2) {
if (val <= 0 | mod(val,1)!= 0) {
errmsg = name + " can only be a positive integer"
_error(3300, errmsg)
}
}
else {
if (val < 0 | mod(val,1)!= 0) {
errmsg = name + " can only be zero or a positive integer"
_error(3300, errmsg)
}
}
}
void clicks::is_pr(real scalar val, string scalar name)
{
string scalar errmsg
if (val > 1 | val < 0) {
errmsg = name + " must remain between 0 and 1"
_error(3300, errmsg)
}
}
transmorphic clicks::Nagents(| real scalar dim)
{
if (args() == 1) {
is_posint(dim, "dim")
Nagents = dim
}
else {
return(Nagents)
}
}
transmorphic clicks::newscycle(| real scalar dim)
{
if (args() == 1) {
is_posint(dim, "dim")
newscycle = dim
}
else {
return(newscycle)
}
}
transmorphic clicks::pr_good(| real scalar dim)
{
if (args() == 1) {
is_pr(dim, "dim")
pr_good = dim
}
else {
return(pr_good)
}
}
transmorphic clicks::tdim(| real scalar dim)
{
if (args() == 1) {
is_posint(dim, "dim")
tdim = dim
}
else {
return(tdim)
}
}
transmorphic clicks::baseclicks(| real scalar dim)
{
if (args() == 1) {
is_posint(dim, "dim")
baseclicks = dim
}
else {
return(baseclicks)
}
}
transmorphic clicks::botclicks(| real scalar dim)
{
if (args() == 1) {
is_posint(dim, "dim", "zero_ok")
botclicks = dim
}
else {
return(botclicks)
}
}
void clicks::setup()
{
if (tdim==.) tdim = 10
if (pr_good==.) pr_good = .5
if (newscycle==.) newscycle = 20
if (Nagents==.) Nagents = 1000
if (baseclicks==.) baseclicks = 2*Nagents/newscycle
if (botclicks==.) botclicks = 0
agents.reinit("real", 2)
items.reinit("real", 3)
pr.reinit("real", 1)
}
void clicks::populate_news()
{
real scalar i
real rowvector key
for(i=1; i <= newscycle ; i++) {
key = i, 0, `content'
items.put(key,runiform(1,1)<pr_good)
key = i, 0, `clicks'
items.put(key,i*ceil(runiform(1,1)*baseclicks))
}
}
void clicks::populate_agents()
{
real scalar i
real rowvector key
for(i=1; i<=Nagents ; i++) {
key = i,0
agents.put(key,ceil(runiform(1,1)*10))
}
}
void clicks::run()
{
real scalar t , i
setup()
populate_news()
populate_agents()
for(t=1; t <= tdim; t++) {
news_step(t)
pr(t)
for(i=1; i <= Nagents; i++) {
agent_step(i,t)
}
}
}
end
mata
foo = clicks()
foo.tdim(200)
foo.pr_good(.5)
foo.botclicks(50)
foo.run()
foo.mood(("m1","m2","m3","m4","m5","m6","m7","m8","m9","m10"))
end
gen t = _n-1
sort t
forvalues i = 2/10 {
local j = `i' - 1
replace m`i' = m`i' + m`j'
}
twoway area m10 t, color("165 0 38") || ///
area m9 t, color("215 48 39") || ///
area m8 t, color("244 109 67") || ///
area m7 t, color("253 174 97") || ///
area m6 t, color("254 224 144") || ///
area m5 t, color("224 243 248") || ///
area m4 t, color("171 217 233") || ///
area m3 t, color("116 173 209") || ///
area m2 t, color(" 69 117 180") || ///
area m1 t, color(" 49 54 149") ///
plotregion(margin(zero)) ///
name(mood, replace) ///
legend(order(1 "10 good" ///
2 "9" ///
3 "8" ///
4 "7" ///
5 "6" ///
6 "5" ///
7 "4" ///
8 "3" ///
9 "2" ///
10 "1 bad" ))
drop _all
mata foo.items(("id", "t", "what", "val"))
reshape wide val, i(id t) j(what)
gen clicks_good = val2*val1
gen clicks_bad = val2*(1-val1)
collapse (sum) clicks_good clicks_bad, by(t)
twoway line clicks_good clicks_bad t, ///
lcolor("165 0 38" " 49 54 149") ///
lpattern(solid solid) ///
legend(order(1 "good" 2 "bad")) ///
ytitle(total number of clicks)
-------------------------------------------------------------------------------
<< index >>
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Conclusion
-------------------------------------------------------------------------------
Conclusion
Implementing and agent based model in Mata is doable
It is a non-trivial programming exercise, but that is to some extend true
for all platforms
Creating a GUI tends to be easier in the other platforms, though that
does not worry me that much. Once the model has been programmed, the
interaction with that model is intuitive enough.
Creating simple tables and graphs also tends to be easier in the other
platforms. I find that more serious.
-------------------------------------------------------------------------------
<< index
-------------------------------------------------------------------------------
